#include "File/YamlLoader.h"
#include "Functions/TerminalIO.h"
#include "Riostream.h"
#include "TROOT.h"
#include "TRolke.h"
#include "TSystem.h"
using namespace std;
using namespace BOSS_Afterburner::TerminalIO;

/// @addtogroup Afterburner_scripts
/// @{

/// Based on the [`Rolke` example](https://nbviewer.jupyter.org/url/root.cern.ch/doc/master/notebooks/Rolke.C.nbconvert.ipynb) that comes with ROOT, but adapted with `YamlLoader`.
/// @see [`TRolke` class description](https://root.cern.ch/doc/master/classTRolke.html).
/// @author   Remco de Boer 雷穆克 (r.e.deboer@students.uu.nl or remco.de.boer@ihep.ac.cn)
/// @date     May 9th, 2019
void Rolke(const char* config_file)
{
  YamlLoader conf(config_file);

  // Number of observed events in the experiment
  Int_t x;

  // Number of observed events in background region
  Int_t y;
  // Either: the ratio between signal and background region in case background is observed.
  // Or: the ratio between observed and simulated live-time in case background is determined from MC.
  Double_t tau;

  // Estimate of the background
  Double_t bm;
  // Standard deviation of the background estimate
  Double_t sdb;
  // Estimate of the efficiency
  Double_t em;
  // Standard deviation of the efficiency estimate
  Double_t sde;

  // Number of MC events generated
  Int_t m;
  // Number of MC events observed
  Int_t z;

  // True efficiency (considered known)
  Double_t e;
  // Background expectation value (considered known)
  Double_t b;
  // Confidence Level
  Double_t alpha;

  // Make TRolke objects
  TRolke   tr; //
  Double_t ul; // Upper limit
  Double_t ll; // Lower limit

  /// Model 1 assumes:
  /// - Poisson uncertainty in the background estimate
  /// - Binomial uncertainty in the efficiency estimate
  const auto& model1 = conf["Model 1"];
  if(model1)
  {
    cout << endl << " ======================================================== " << endl;
    x     = model1["Number of observed events in the experiment"].as<Int_t>();
    y     = model1["Number of observed events in background region"].as<Int_t>();
    tau   = model1["Ratio between signal and background region"].as<Double_t>();
    m     = model1["Number of MC events generated"].as<Int_t>();
    z     = model1["Number of MC events observed"].as<Int_t>();
    alpha = model1["Confidence Level"].as<Double_t>();
    tr.SetCL(alpha);
    tr.SetPoissonBkgBinomEff(x, y, z, tau, m);
    tr.GetLimits(ll, ul);
    cout << "For model 1: Poisson / Binomial" << endl;
    cout << "the Profile Likelihood interval is :" << endl;
    cout << "[" << ll << "," << ul << "]" << endl;
  }

  /// Model 2 assumes:
  /// - Poisson uncertainty in the background estimate
  /// - Gaussian  uncertainty in the efficiency estimate
  const auto& model2 = conf["Model 2"];
  if(model2)
  {
    cout << endl << " ======================================================== " << endl;
    x     = model2["Number of observed events in the experiment"].as<Int_t>();
    y     = model2["Number of observed events in background region"].as<Int_t>();
    tau   = model2["Ratio between signal and background region"].as<Double_t>();
    em    = model2["Estimate of the efficiency"].as<Double_t>();
    sde   = model2["Standard deviation of the efficiency estimate"].as<Double_t>();
    alpha = model2["Confidence Level"].as<Double_t>();
    tr.SetCL(alpha);
    tr.SetPoissonBkgGaussEff(x, y, em, tau, sde);
    tr.GetLimits(ll, ul);
    cout << "For model 2 : Poisson / Gaussian" << endl;
    cout << "the Profile Likelihood interval is :" << endl;
    cout << "[" << ll << "," << ul << "]" << endl;
  }

  /// Model 3 assumes:
  /// - Gaussian uncertainty in the background estimate
  /// - Gaussian  uncertainty in the efficiency estimate
  const auto& model3 = conf["Model 3"];
  if(model3)
  {
    cout << endl << " ======================================================== " << endl;
    x     = model3["Number of observed events in the experiment"].as<Int_t>();
    bm    = model3["Estimate of the background"].as<Double_t>();
    sdb   = model3["Standard deviation of the background estimate"].as<Double_t>();
    em    = model3["Estimate of the efficiency"].as<Double_t>();
    sde   = model3["Standard deviation of the efficiency estimate"].as<Double_t>();
    alpha = model3["Confidence Level"].as<Double_t>();
    tr.SetCL(alpha);
    tr.SetGaussBkgGaussEff(x, bm, em, sde, sdb);
    tr.GetLimits(ll, ul);
    cout << "For model 3 : Gaussian / Gaussian" << endl;
    cout << "the Profile Likelihood interval is :" << endl;
    cout << "[" << ll << "," << ul << "]" << endl;
    cout << "***************************************" << endl;
    cout << "* some more example's for gauss/gauss *" << endl;
    cout << "*                                     *" << endl;
    Double_t slow, shigh;
    tr.GetSensitivity(slow, shigh);
    cout << "sensitivity:" << endl;
    cout << "[" << slow << "," << shigh << "]" << endl;
    int outx;
    tr.GetLimitsQuantile(slow, shigh, outx, 0.5);
    cout << "median limit:" << endl;
    cout << "[" << slow << "," << shigh << "] @ x =" << outx << endl;
    tr.GetLimitsML(slow, shigh, outx);
    cout << "ML limit:" << endl;
    cout << "[" << slow << "," << shigh << "] @ x =" << outx << endl;
    tr.GetSensitivity(slow, shigh);
    cout << "sensitivity:" << endl;
    cout << "[" << slow << "," << shigh << "]" << endl;
    tr.GetLimits(ll, ul);
    cout << "the Profile Likelihood interval is :" << endl;
    cout << "[" << ll << "," << ul << "]" << endl;
    Int_t ncrt;
    tr.GetCriticalNumber(ncrt);
    cout << "critical number: " << ncrt << endl;
    tr.SetCLSigmas(5);
    tr.GetCriticalNumber(ncrt);
    cout << "critical number for 5 sigma: " << ncrt << endl;
    cout << "***************************************" << endl;
  }

  /// Model 4 assumes:
  /// - Poisson uncertainty in the background estimate
  /// - known efficiency
  const auto& model4 = conf["Model 4"];
  if(model4)
  {
    cout << endl << " ======================================================== " << endl;
    x     = model4["Number of observed events in the experiment"].as<Int_t>();
    y     = model4["Number of observed events in background region"].as<Int_t>();
    tau   = model4["Ratio between signal and background region"].as<Double_t>();
    e     = model4["True efficiency (considered known)"].as<Double_t>();
    alpha = model4["Confidence Level"].as<Double_t>();
    tr.SetCL(alpha);
    tr.SetPoissonBkgKnownEff(x, y, tau, e);
    tr.GetLimits(ll, ul);
    cout << "For model 4 : Poissonian / Known" << endl;
    cout << "the Profile Likelihood interval is :" << endl;
    cout << "[" << ll << "," << ul << "]" << endl;
  }

  /// Model 5 assumes:
  /// - Gaussian uncertainty in the background estimate
  /// - Known efficiency
  const auto& model5 = conf["Model 5"];
  if(model5)
  {
    cout << endl << " ======================================================== " << endl;
    bm    = model5["Estimate of the background"].as<Double_t>();
    x     = model5["Number of observed events in the experiment"].as<Int_t>();
    e     = model5["True efficiency (considered known)"].as<Double_t>();
    sdb   = model5["Standard deviation of background estimate"].as<Double_t>();
    alpha = model5["Confidence Level"].as<Double_t>();
    tr.SetCL(alpha);
    tr.SetGaussBkgKnownEff(x, bm, sdb, e);
    tr.GetLimits(ll, ul);
    cout << "For model 5 : Gaussian / Known" << endl;
    cout << "the Profile Likelihood interval is :" << endl;
    cout << "[" << ll << "," << ul << "]" << endl;
  }

  /// Model 6 assumes:
  /// - Background expectation value (considered known)
  /// - Binomial uncertainty in the efficiency estimate
  const auto& model6 = conf["Model 6"];
  if(model6)
  {
    cout << endl << " ======================================================== " << endl;
    b     = model6["Background expectation value (considered known)"].as<Double_t>();
    x     = model6["Number of observed events in the experiment"].as<Int_t>();
    m     = model6["Number of MC events generated"].as<Int_t>();
    z     = model6["Number of MC events observed"].as<Int_t>();
    alpha = model6["Confidence Level"].as<Double_t>();
    tr.SetCL(alpha);
    tr.SetKnownBkgBinomEff(x, z, m, b);
    tr.GetLimits(ll, ul);
    cout << "For model 6 : Known / Binomial" << endl;
    cout << "the Profile Likelihood interval is :" << endl;
    cout << "[" << ll << "," << ul << "]" << endl;
  }

  /// Model 7 assumes:
  /// - Known Background
  /// - Gaussian  uncertainty in the efficiency estimate
  const auto& model7 = conf["Model 7"];
  if(model7)
  {
    cout << endl << " ======================================================== " << endl;
    x     = model7["Number of observed events in the experiment"].as<Int_t>();
    em    = model7["Estimate of the efficiency"].as<Double_t>();
    sde   = model7["Standard deviation of the efficiency estimate"].as<Double_t>();
    b     = model7["Background expectation value (considered known)"].as<Double_t>();
    alpha = model7["Confidence Level"].as<Double_t>();
    y     = 1;
    tr.SetCL(alpha);
    tr.SetKnownBkgGaussEff(x, em, sde, b);
    tr.GetLimits(ll, ul);
    cout << "For model 7 : Known / Gaussian " << endl;
    cout << "the Profile Likelihood interval is :" << endl;
    cout << "[" << ll << "," << ul << "]" << endl;
  }

  /// # Example of bounded and unbounded likelihood
  /// Example for Model 1
  const auto& bounding = conf["Bounding"];
  if(bounding)
  {
    bm    = bounding["Estimate of the background"].as<Double_t>();
    tau   = bounding["Ratio between signal and background region"].as<Double_t>();
    m     = bounding["Number of MC events generated"].as<Int_t>();
    z     = bounding["Number of MC events observed"].as<Int_t>();
    y     = bounding["Number of observed events in background region"].as<Int_t>();
    x     = bounding["Number of observed events in the experiment"].as<Int_t>();
    alpha = bounding["Confidence Level"].as<Double_t>();
    tr.SetCL(alpha);
    tr.SetPoissonBkgBinomEff(x, y, z, tau, m);
    tr.SetBounding(true); // bounded
    tr.GetLimits(ll, ul);
    cout << "Example of the effect of bounded vs unbounded, For model 1" << endl;
    cout << "the BOUNDED Profile Likelihood interval is :" << endl;
    cout << "[" << ll << "," << ul << "]" << endl;
    tr.SetBounding(false); // unbounded
    tr.GetLimits(ll, ul);
    cout << "the UNBOUNDED Profile Likelihood interval is :" << endl;
    cout << "[" << ll << "," << ul << "]" << endl;
  }
}
/// @}

int main(int argc, char* argv[])
{
  CheckMainArguments(argc, argv, 1);
  Rolke(argv[1]);
  return 0;
}